Stress x EE MA
01/07/2021
Loading packages
knitr::opts_chunk$set(echo = TRUE)
pacman::p_load(tidyverse,
here,
metafor,
emmeans,
orchaRd,
MCMCglmm,
corpcor)Loading data and functions
dat <- read_csv(here("Data","Full_data.csv"))
# Load custom function to extract data
source(here("R/functions.R")) Data organisation
Getting effect sizes from function, ‘flipping’ effect sizes so that all effect sizes are higher values = individuals do better and learning/memory, and shifting negative values to possitive as lnRR cannot use negative values
#removing negative values
#dat <- dat %>%
# mutate(CC_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
# CS_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
# EC_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
# ES_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
# CC_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
# CS_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
# EC_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
# ES_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
#CC_mean = case_when(ES_ID == 38 ~ CC_mean+80),
#CS_mean = case_when(ES_ID == 38 ~ CC_mean+80),
#EC_mean = case_when(ES_ID == 38 ~ CC_mean+80),
#ES_mean = case_when(ES_ID == 38 ~ CC_mean+80),
#CC_mean = case_when(ES_ID == 39 ~ CC_mean+33.86),
#CS_mean = case_when(ES_ID == 39 ~ CC_mean+33.86),
#EC_mean = case_when(ES_ID == 39 ~ CC_mean+33.86),
#ES_mean = case_when(ES_ID == 39 ~ CC_mean+33.86))
#dat$CC_mean <- ifelse(dat$ES_ID == 36, dat$CC_mean+29.4, dat$CC_mean)
#Getting effect sizes
effect_size <- effect_set(CC_n = "CC_n", CC_mean = "CC_mean", CC_SD = "CC_SD",
EC_n = "EC_n", EC_mean = "EC_mean" , EC_SD ="EC_SD",
CS_n = "CS_n", CS_mean = "CS_mean", CS_SD = "CS_SD",
ES_n = "ES_n", ES_mean = "ES_mean", ES_SD = "ES_SD",
data = dat)
#'pure' effect sizes
effect_size2 <- effect_set2(CC_n = "CC_n", CC_mean = "CC_mean", CC_SD = "CC_SD",
EC_n = "EC_n", EC_mean = "EC_mean" , EC_SD ="EC_SD",
CS_n = "CS_n", CS_mean = "CS_mean", CS_SD = "CS_SD",
ES_n = "ES_n", ES_mean = "ES_mean", ES_SD = "ES_SD",
data = dat)
#Removing missing effect sizes
dim(dat)
full_info <- which(complete.cases(effect_size) == TRUE)
dat_effect <- cbind(dat, effect_size, effect_size2)
dat <- dat_effect[full_info, ]
names(dat_effect)
# TODO - we need to sort this out a bit more later
#There are no longer any missing effect sizes but is missing info in moderators: do we still need to remove these?
dat <- dat_effect[full_info, ]
dim(dat_effect)
dimentions <- dim(dat) # 7 less
#Flipping 'lower is better' effect sizes
#flipping lnRR for values where higher = worse
dat$lnRR_Ea <- ifelse(dat$Response_direction == 2, dat$lnRR_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E)) # currently NAswhich causes error
dat$lnRR_Sa <- ifelse(dat$Response_direction == 2, dat$lnRR_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S)) # currently NAswhich causes error
dat$lnRR_ESa <- ifelse(dat$Response_direction == 2, dat$lnRR_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES)) # currently NAswhich causes error
#flipping 'pure effect sizes'
dat$lnRR_E2a <- ifelse(dat$Response_direction == 2, dat$lnRR_E2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E2)) # currently NAswhich causes error
dat$lnRR_S2a <- ifelse(dat$Response_direction == 2, dat$lnRR_S2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S2)) # currently NAswhich causes error
dat$lnRR_ES2a <- ifelse(dat$Response_direction == 2, dat$lnRR_ES2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES2)) # currently NAswhich causes error
dat$lnRR_E3a <- ifelse(dat$Response_direction == 2, dat$lnRR_E3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E3)) # currently NAswhich causes error
dat$lnRR_S3a <- ifelse(dat$Response_direction == 2, dat$lnRR_S3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S3)) # currently NAswhich causes error
#flipping SMD
dat$SMD_Ea <- ifelse(dat$Response_direction == 2, dat$SMD_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_E)) # currently NAswhich causes error
dat$SMD_Sa <- ifelse(dat$Response_direction == 2, dat$SMD_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_S)) # currently NAswhich causes error
dat$SMD_ESa <- ifelse(dat$Response_direction == 2, dat$SMD_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_ES))
#rounding down integer sample sizes
dat$CC_n <- floor(dat$CC_n)
dat$EC_n <- floor(dat$EC_n)
dat$CS_n <- floor(dat$CS_n)
dat$ES_n <- floor(dat$CS_n)Assigning human readable words to moderators
dat <- dat %>% mutate(Type_learning = case_when(Type_learning == 1 ~ "Habituation",
Type_learning == 2 ~ "Conditioning",
Type_learning == 3 ~ "Recognition",
Type_learning == 4 ~ "Unclear"),
Learning_vs_memory = case_when(Learning_vs_memory == 1 ~ "Learning",
Learning_vs_memory == 2 ~ "Memory",
Learning_vs_memory == 1 ~ "Unclear"),
Appetitive_vs_aversive = case_when(Appetitive_vs_aversive == 1 ~"Appetitive",
Appetitive_vs_aversive == 2 ~ "Aversive",
Appetitive_vs_aversive == 3 ~ "Not applicable",
Appetitive_vs_aversive == 4 ~ "Unclear"),
Type_stress_exposure = case_when(Type_stress_exposure == 1 ~ "Density",
Type_stress_exposure == 2 ~ "Scent",
Type_stress_exposure == 3 ~ "Shock",
Type_stress_exposure == 4 ~ "Exertion",
Type_stress_exposure == 5 ~ "Restraint",
Type_stress_exposure == 6 ~ "MS",
Type_stress_exposure == 7 ~ "Circadian rhythm",
Type_stress_exposure == 8 ~ "Noise",
Type_stress_exposure == 9 ~ "Other",
Type_stress_exposure == 10 ~ "Combination",
Type_stress_exposure == 11 ~ "unclear"),
Age_stress_exposure = case_when(Age_stress_exposure == 1 ~ "Prenatal",
Age_stress_exposure == 2 ~ "Juvenile",
Age_stress_exposure == 3 ~ "Adult",
Age_stress_exposure == 4 ~ "Unclear"),
Stress_duration = case_when(Stress_duration == 1 ~ "Acute",
Stress_duration == 2 ~ "Chronic",
Stress_duration == 3 ~ "Intermittent",
Stress_duration == 4 ~ "Unclear"),
EE_social = case_when(EE_social == 1 ~ "Social",
EE_social== 2 ~ "Non-social",
EE_social == 3 ~ "Unclear"),
EE_exercise = case_when(EE_exercise == 1 ~ "Exercise",
EE_exercise == 2 ~ "No exercise"),
Age_EE_exposure = case_when(Age_EE_exposure == 1 ~ "Prenatal",
Age_EE_exposure == 2 ~ "Juvenile",
Age_EE_exposure == 3 ~ "Adult",
Age_EE_exposure == 4 ~ "Unclear"))Making Variance VCV
Used to control for trait relatedness within individuals
##adding data to empty VCV matrix
VCV_E <- make_VCV_matrix(data = dat, V = "lnRRV_E", cluster = "Study_ID", obs = "ES_ID")
VCV_S <- make_VCV_matrix(data = dat, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")
VCV_ES <- make_VCV_matrix(data = dat, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
is.positive.definite(VCV_E)
isSymmetric(VCV_E)
is.positive.definite(VCV_S)
isSymmetric(VCV_S)
is.positive.definite(VCV_ES)
isSymmetric(VCV_ES)
VCV_Ea <- make_VCV_matrix(data = dat, V = "SMDV_E", cluster = "Study_ID", obs = "ES_ID")
VCV_Sa <- make_VCV_matrix(data = dat, V = "SMDV_S", cluster = "Study_ID", obs = "ES_ID")
VCV_ESa <- make_VCV_matrix(data = dat, V = "SMDV_ES", cluster = "Study_ID", obs = "ES_ID")
is.positive.definite(VCV_Ea)
isSymmetric(VCV_Ea)
is.positive.definite(VCV_Sa)
isSymmetric(VCV_Sa)
is.positive.definite(VCV_ESa)
isSymmetric(VCV_ESa)Examining the group of large/positive effect sizes
Note that I checked the top four most positive interaction effect sizes and the data looks correct. Highest ES is study 27, effect size 88, where Maternally separated EE individuals freeze significantly more when given a cue from inhibitory avoidance. Im not sure why the second highest one (ES 16) is so high as the control indivs do the best. EE indivs do best in the 3rd highest ES (ES 10) so not sure about this one either? The fourth highest (ES 12) seems correct
largeES_S <- dat %>%
slice_max(lnRR_S, n = 10)
largeES_E <- dat %>%
slice_max(lnRR_E, n = 10)
largeES_ES <-dat %>%
slice_max(lnRR_ES, n = 10)
#most of the largest effect sizes come from different papersModelling with lnRR
Things todo: - Incorporate strain as random effect - Consider including VCV - Check outlier - Remove moderator levels with k < 5
Stress
Intercept model
Learning and memory are significantly reduced due to stress. High heterogeneity
mod_S0 <- rma.mv(yi = lnRR_Sa, V = VCV_S, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.5945 21.1890 29.1890 39.2325 29.6542
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0080 0.0892 30 no Study_ID
## sigma^2.2 0.0358 0.1891 92 no ES_ID
## sigma^2.3 0.0001 0.0086 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 863.4222, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0903 0.0331 -2.7240 91 0.0077 -0.1561 -0.0244 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_S0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 0.928503647 0.168668323 0.758268848 0.001566476funnel(mod_S0)
orchard_plot(mod_S0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Publication Bias
… I was going to run this using the same code as last meta-analysis (MCMCglmm) but thought there might be an easier way?
#need to remove NAs
dat_na_S <- dat %>%
drop_na(c(lnRR_Sa, lnRRV_S))
dat_na_E <- dat %>%
drop_na(c(lnRR_Ea, lnRRV_E))
dat_na_ES <- dat %>%
drop_na(c(lnRR_ESa, lnRRV_ES))
#Stress
#TODO check with Shinichi if there is an easier way to do PB (eggers regression, funnel asymmetry)Meta-regression
Type of learning
The type of learning/memory response
dat$Type_learning<-as.factor(dat$Type_learning)
dat1 <- filter(dat, Type_learning %in% c("Recognition", "Habituation", "Conditioning"))
VCV_S1 <- make_VCV_matrix(data = dat1, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")
mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_learning-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_S1)##
## Multivariate Meta-Analysis Model (k = 91; method: REML)
##
## logLik Deviance AIC BIC AICc
## -4.8116 9.6232 21.6232 36.4872 22.6602
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0196 0.1400 30 no Study_ID
## sigma^2.2 0.0192 0.1385 91 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 498.8005, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 5.2846, p-val = 0.0021
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_learningConditioning -0.1110 0.0387 -2.8663 88 0.0052 -0.1880
## Type_learningHabituation -0.2548 0.0816 -3.1212 88 0.0024 -0.4170
## Type_learningRecognition 0.0096 0.0571 0.1686 88 0.8665 -0.1038
## ci.ub
## Type_learningConditioning -0.0340 **
## Type_learningHabituation -0.0926 **
## Type_learningRecognition 0.1230
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S1) ## R2_marginal R2_coditional
## 0.08953573 1.00000000# Orchard plot
orchard_plot(mod_S1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Learning vs Memory
Is the assay broadly measuring learning or memory?
mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S2) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.8490 21.6980 31.6980 43.7922 32.4772
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0184 0.1357 30 no Study_ID
## sigma^2.2 0.0279 0.1669 85 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 623.2963, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 2.7188, p-val = 0.0718
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning -0.0342 0.0509 -0.6715 83 0.5037 -0.1354
## Learning_vs_memoryMemory -0.0922 0.0400 -2.3073 83 0.0235 -0.1718
## ci.ub
## Learning_vs_memoryLearning 0.0670
## Learning_vs_memoryMemory -0.0127 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S2) ## R2_marginal R2_coditional
## 0.01603098 1.00000000# Orchard plot
orchard_plot(mod_S2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Appetitive vs aversive
The type of cue used
dat2 <- filter(dat, Appetitive_vs_aversive %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_S2 <- make_VCV_matrix(data = dat2, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")
mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~ Appetitive_vs_aversive-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_S3)##
## Multivariate Meta-Analysis Model (k = 91; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.5287 13.0575 25.0575 39.9215 26.0945
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0145 0.1203 30 no Study_ID
## sigma^2.2 0.0256 0.1601 91 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 506.6327, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 3.8651, p-val = 0.0120
##
## Model Results:
##
## estimate se tval df pval
## Appetitive_vs_aversiveAppetitive -0.1781 0.0756 -2.3563 88 0.0207
## Appetitive_vs_aversiveAversive -0.1065 0.0429 -2.4801 88 0.0150
## Appetitive_vs_aversiveNot applicable -0.0320 0.0508 -0.6306 88 0.5300
## ci.lb ci.ub
## Appetitive_vs_aversiveAppetitive -0.3283 -0.0279 *
## Appetitive_vs_aversiveAversive -0.1919 -0.0212 *
## Appetitive_vs_aversiveNot applicable -0.1329 0.0689
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S3) ## R2_marginal R2_coditional
## 0.05569185 1.00000000# Orchard plot
orchard_plot(mod_S3, mod = "Appetitive_vs_aversive", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Type of stress
The type of stress manipulation
dat3 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"))
VCV_S3 <- make_VCV_matrix(data = dat3, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")
mod_S4 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_S4) ##
## Multivariate Meta-Analysis Model (k = 87; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.2101 18.4202 32.4202 49.3521 33.9135
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0080 0.0896 26 no Study_ID
## sigma^2.2 0.0373 0.1932 87 no ES_ID
## sigma^2.3 0.0000 0.0000 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 764.1363, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 83) = 2.7813, p-val = 0.0319
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination -0.0809 0.0851 -0.9511 83 0.3443 -0.2502
## Type_stress_exposureMS -0.0507 0.0508 -0.9983 83 0.3210 -0.1518
## Type_stress_exposureNoise -0.0739 0.0905 -0.8164 83 0.4166 -0.2539
## Type_stress_exposureRestraint -0.2026 0.0692 -2.9253 83 0.0044 -0.3403
## ci.ub
## Type_stress_exposureCombination 0.0883
## Type_stress_exposureMS 0.0503
## Type_stress_exposureNoise 0.1061
## Type_stress_exposureRestraint -0.0648 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S4)## R2_marginal R2_coditional
## 0.0690348 1.0000000# Orchard plot
orchard_plot(mod_S4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Age of stress
The age at which the individuals were exposed to the stressor. Note there are a lot of ‘unkown’ age as authors only report PND which needs to be researched. I’m wondering if this also needs an ‘adolescence’ category as this seesm to be popular in rodent research
#TODO these ages need to be recategoried and analysis run again
mod_S5 <-rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S5) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.2914 18.5827 32.5827 49.9241 33.9827
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0014 0.0372 30 no Study_ID
## sigma^2.2 0.0381 0.1953 92 no ES_ID
## sigma^2.3 0.0032 0.0565 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 800.8761, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.9490, p-val = 0.1094
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_stress_exposureAdult -0.1467 0.0857 -1.7123 88 0.0904 -0.3171
## Age_stress_exposureJuvenile -0.0012 0.0543 -0.0214 88 0.9829 -0.1090
## Age_stress_exposurePrenatal -0.1276 0.0873 -1.4617 88 0.1474 -0.3011
## Age_stress_exposureUnclear -0.1588 0.0874 -1.8172 88 0.0726 -0.3325
## ci.ub
## Age_stress_exposureAdult 0.0236 .
## Age_stress_exposureJuvenile 0.1067
## Age_stress_exposurePrenatal 0.0459
## Age_stress_exposureUnclear 0.0149 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S5) ## R2_marginal R2_coditional
## 0.1009404 0.9327030# Orchard plot
orchard_plot(mod_S5, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Acute vs chronic stress
How long was the stress applied for (chronic = every day for 7 days or more)? This has the highest marginal R2
dat4 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- make_VCV_matrix(data = dat3, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")
mod_S6 <-rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat4)
summary(mod_S6) ##
## Multivariate Meta-Analysis Model (k = 87; method: REML)
##
## logLik Deviance AIC BIC AICc
## -20.6635 41.3269 51.3269 63.5402 52.0864
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0085 0.0922 27 no Study_ID
## sigma^2.2 0.0700 0.2645 87 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 85) = 2809.8421, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 85) = 6.9113, p-val = 0.0017
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute 0.1099 0.0873 1.2581 85 0.2118 -0.0638 0.2835
## Stress_durationChronic -0.1569 0.0453 -3.4622 85 0.0008 -0.2470 -0.0668
##
## Stress_durationAcute
## Stress_durationChronic ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S6) ## R2_marginal R2_coditional
## 0.135429 1.000000# Orchard plot
orchard_plot(mod_S6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Environmental enrichment
Intercept model
Learning and memory are significantly improved when housed with environmnetal enrichment
mod_E0 <- rma.mv(yi = lnRR_Ea, V = VCV_E, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.9611 27.9223 35.9223 45.9657 36.3874
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0049 0.0698 30 no Study_ID
## sigma^2.2 0.0428 0.2070 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 914.4170, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1891 0.0321 5.8832 91 <.0001 0.1253 0.2530 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_E0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.339783e-01 9.534725e-02 8.386310e-01 1.741556e-09funnel(mod_E0)
#trying orchard plot
orchard_plot(mod_E0, mod = "Int", xlab = "lnRR", alpha=0.4) + # Orchard plot
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = "darkorange")+ # change colours
scale_fill_manual(values="darkorange")+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Meta-regression
Type of learning
The type of learning/memory response
VCV_E1 <- make_VCV_matrix(data = dat1, V = "lnRRV_E", cluster = "Study_ID", obs = "ES_ID")
mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_learning-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_E1)##
## Multivariate Meta-Analysis Model (k = 91; method: REML)
##
## logLik Deviance AIC BIC AICc
## -11.7462 23.4925 35.4925 50.3565 36.5295
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0145 0.1205 30 no Study_ID
## sigma^2.2 0.0346 0.1859 91 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 823.6996, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 12.2253, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_learningConditioning 0.2415 0.0403 5.9907 88 <.0001 0.1614
## Type_learningHabituation 0.1760 0.0948 1.8565 88 0.0667 -0.0124
## Type_learningRecognition 0.0534 0.0654 0.8167 88 0.4163 -0.0765
## ci.ub
## Type_learningConditioning 0.3217 ***
## Type_learningHabituation 0.3643 .
## Type_learningRecognition 0.1833
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E1) ## R2_marginal R2_coditional
## 0.08695344 1.00000000# Orchard plot
orchard_plot(mod_E1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Learning vs Memory
Is the assay broadly measuring learning or memory?
mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E2) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -8.9742 17.9484 27.9484 40.0426 28.7276
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0074 0.0860 30 no Study_ID
## sigma^2.2 0.0291 0.1707 85 no ES_ID
## sigma^2.3 0.0006 0.0255 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 580.2953, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 13.4308, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.1991 0.0505 3.9454 83 0.0002 0.0987
## Learning_vs_memoryMemory 0.1837 0.0392 4.6870 83 <.0001 0.1057
## ci.ub
## Learning_vs_memoryLearning 0.2994 ***
## Learning_vs_memoryMemory 0.2616 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E2) ## R2_marginal R2_coditional
## 0.001428381 0.982584068# Orchard plot
orchard_plot(mod_E2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Appetitive vs aversive
The type of cue used
VCV_E2 <- make_VCV_matrix(data = dat2, V = "lnRRV_E", cluster = "Study_ID", obs = "ES_ID")
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~ Appetitive_vs_aversive-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_E3)##
## Multivariate Meta-Analysis Model (k = 91; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.0470 18.0940 30.0940 44.9581 31.1311
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0132 0.1148 30 no Study_ID
## sigma^2.2 0.0308 0.1755 91 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 638.2083, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 14.8231, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval
## Appetitive_vs_aversiveAppetitive 0.1895 0.0762 2.4855 88 0.0148
## Appetitive_vs_aversiveAversive 0.2719 0.0440 6.1808 88 <.0001
## Appetitive_vs_aversiveNot applicable 0.0623 0.0528 1.1792 88 0.2415
## ci.lb ci.ub
## Appetitive_vs_aversiveAppetitive 0.0380 0.3410 *
## Appetitive_vs_aversiveAversive 0.1845 0.3594 ***
## Appetitive_vs_aversiveNot applicable -0.0427 0.1673
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E3) ## R2_marginal R2_coditional
## 0.1467484 1.0000000# Orchard plot
orchard_plot(mod_E3, mod = "Appetitive_vs_aversive", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Exercise enrichment
Does the form of enrichment include exercise through a running wheel or treadmill?
mod_E5<- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.1841 28.3682 38.3682 50.8673 39.0825
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0055 0.0744 30 no Study_ID
## sigma^2.2 0.0431 0.2077 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 889.4288, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.8895, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1868 0.0396 4.7137 90 <.0001 0.1081 0.2656
## EE_exerciseNo exercise 0.1953 0.0574 3.3999 90 0.0010 0.0812 0.3094
##
## EE_exerciseExercise ***
## EE_exerciseNo exercise **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E5) ## R2_marginal R2_coditional
## 0.0003377362 0.9999999972# Orchard plot
orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Age of enrichment
The age at which the individuals were exposed to environmental enrichment. This needs to be redone with new age categories
mod_E6 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E6) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -12.7731 25.5462 37.5462 52.4781 38.5706
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0049 0.0699 30 no Study_ID
## sigma^2.2 0.0437 0.2090 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 882.6717, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.5329, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Age_EE_exposureAdult 0.1838 0.0648 2.8351 89 0.0057 0.0550 0.3126
## Age_EE_exposureJuvenile 0.0082 0.1072 0.0767 89 0.9390 -0.2047 0.2212
## Age_EE_exposureUnclear 0.2164 0.0398 5.4365 89 <.0001 0.1373 0.2955
##
## Age_EE_exposureAdult **
## Age_EE_exposureJuvenile
## Age_EE_exposureUnclear ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E6) ## R2_marginal R2_coditional
## 0.07138796 1.00000000# Orchard plot
orchard_plot(mod_E6, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Interaction of stress and EE
Intercept
Enriched and stress animals are better at learning and memory. TODO: It looks like there is a large but low precision outlier. Should potentially remove?
mod_ES0 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -45.8577 91.7153 99.7153 109.7588 100.1804
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0347 0.1862 30 no Study_ID
## sigma^2.2 0.0269 0.1640 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 322.9690, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1288 0.0496 2.5996 91 0.0109 0.0304 0.2273 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_ES0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.202982e-01 4.617725e-01 3.585257e-01 1.352473e-08funnel(mod_ES0)
orchard_plot(mod_ES0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
#running the model again with the extreme outlier removed
dat_outliers <- dat %>%
filter(ES_ID != 88)
VCV_ES_outliers <- make_VCV_matrix(data = dat_outliers, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ESoutlier <- rma.mv(yi = lnRR_ESa, V = VCV_ES_outliers, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_outliers)
summary(mod_ES0)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -45.8577 91.7153 99.7153 109.7588 100.1804
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0347 0.1862 30 no Study_ID
## sigma^2.2 0.0269 0.1640 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 322.9690, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1288 0.0496 2.5996 91 0.0109 0.0304 0.2273 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1orchard_plot(mod_ESoutlier, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Meta-regression
Type of learning
The type of learning/memory response
VCV_ES1 <- make_VCV_matrix(data = dat1, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_learning-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_ES1)##
## Multivariate Meta-Analysis Model (k = 91; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.8899 85.7798 97.7798 112.6438 98.8169
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0390 0.1974 30 no Study_ID
## sigma^2.2 0.0194 0.1392 91 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 306.3683, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 5.3955, p-val = 0.0019
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_learningConditioning 0.1874 0.0550 3.4087 88 0.0010 0.0782
## Type_learningHabituation 0.1787 0.1069 1.6724 88 0.0980 -0.0336
## Type_learningRecognition -0.0429 0.0740 -0.5800 88 0.5634 -0.1900
## ci.ub
## Type_learningConditioning 0.2967 ***
## Type_learningHabituation 0.3911 .
## Type_learningRecognition 0.1041
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES1) ## R2_marginal R2_coditional
## 0.1063433 1.0000000# Orchard plot
orchard_plot(mod_ES1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Learning vs Memory
Is the assay broadly measuring learning or memory?
mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES2) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -47.7835 95.5669 105.5669 117.6611 106.3461
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0345 0.1857 30 no Study_ID
## sigma^2.2 0.0313 0.1769 85 no ES_ID
## sigma^2.3 0.0000 0.0001 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 313.8032, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 3.7300, p-val = 0.0281
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.1817 0.0694 2.6157 83 0.0106 0.0435
## Learning_vs_memoryMemory 0.1066 0.0540 1.9724 83 0.0519 -0.0009
## ci.ub
## Learning_vs_memoryLearning 0.3198 *
## Learning_vs_memoryMemory 0.2141 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES2) ## R2_marginal R2_coditional
## 0.0187792 0.9999999# Orchard plot
orchard_plot(mod_ES2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Appetitive vs aversive
The type of cue used
VCV_ES2 <- make_VCV_matrix(data = dat2, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~ Appetitive_vs_aversive-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_ES3)##
## Multivariate Meta-Analysis Model (k = 91; method: REML)
##
## logLik Deviance AIC BIC AICc
## -43.0121 86.0243 98.0243 112.8883 99.0613
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0327 0.1809 30 no Study_ID
## sigma^2.2 0.0218 0.1476 91 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 281.9381, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 5.0524, p-val = 0.0028
##
## Model Results:
##
## estimate se tval df pval
## Appetitive_vs_aversiveAppetitive 0.1612 0.1081 1.4909 88 0.1396
## Appetitive_vs_aversiveAversive 0.1986 0.0594 3.3409 88 0.0012
## Appetitive_vs_aversiveNot applicable 0.0126 0.0652 0.1929 88 0.8475
## ci.lb ci.ub
## Appetitive_vs_aversiveAppetitive -0.0537 0.3761
## Appetitive_vs_aversiveAversive 0.0805 0.3168 **
## Appetitive_vs_aversiveNot applicable -0.1169 0.1421
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES3) ## R2_marginal R2_coditional
## 0.09870464 1.00000000# Orchard plot
orchard_plot(mod_ES3, mod = "Appetitive_vs_aversive", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Type of stress
The type of stress manipulation
VCV_ES3 <- make_VCV_matrix(data = dat3, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_ES4)##
## Multivariate Meta-Analysis Model (k = 87; method: REML)
##
## logLik Deviance AIC BIC AICc
## -41.9621 83.9243 97.9243 114.8562 99.4176
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0469 0.2167 26 no Study_ID
## sigma^2.2 0.0274 0.1655 87 no ES_ID
## sigma^2.3 0.0000 0.0000 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 292.2191, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 83) = 1.3135, p-val = 0.2718
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination 0.0995 0.1398 0.7121 83 0.4784 -0.1785
## Type_stress_exposureMS 0.1167 0.0853 1.3684 83 0.1749 -0.0529
## Type_stress_exposureNoise 0.1709 0.1592 1.0735 83 0.2862 -0.1458
## Type_stress_exposureRestraint 0.1480 0.1128 1.3122 83 0.1931 -0.0764
## ci.ub
## Type_stress_exposureCombination 0.3775
## Type_stress_exposureMS 0.2863
## Type_stress_exposureNoise 0.4876
## Type_stress_exposureRestraint 0.3724
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES4)## R2_marginal R2_coditional
## 0.007436653 0.999999995# Orchard plot
orchard_plot(mod_ES4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Age of stress
The age at which the individuals were exposed to the stressor. Note there are a lot of ‘unkown’ age as authors only report PND which needs to be researched. I’m wondering if this also needs an ‘adolescence’ category as this seesm to be popular in rodent research This nees to be redone after age categories are re done
mod_ES5 <-rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES5) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.6838 85.3676 99.3676 116.7089 100.7676
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0238 0.1542 30 no Study_ID
## sigma^2.2 0.0259 0.1610 92 no ES_ID
## sigma^2.3 0.0131 0.1146 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 282.2171, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 2.2788, p-val = 0.0671
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_stress_exposureAdult 0.0739 0.1358 0.5438 88 0.5880 -0.1960
## Age_stress_exposureJuvenile 0.0472 0.0941 0.5008 88 0.6177 -0.1399
## Age_stress_exposurePrenatal 0.3611 0.1492 2.4204 88 0.0176 0.0646
## Age_stress_exposureUnclear 0.2506 0.1345 1.8624 88 0.0659 -0.0168
## ci.ub
## Age_stress_exposureAdult 0.3437
## Age_stress_exposureJuvenile 0.2342
## Age_stress_exposurePrenatal 0.6576 *
## Age_stress_exposureUnclear 0.5179 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES5) ## R2_marginal R2_coditional
## 0.1611826 0.8245990# Orchard plot
orchard_plot(mod_ES5, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Acute vs chronic stress
How long was the stress applied for (chronic = every day for 7 days or more)? This has the highest marginal R2 (currentl nearly 43%) - need to redo without outlier
VCV_ES4 <- make_VCV_matrix(data = dat4, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES6 <-rma.mv(yi = lnRR_ESa, V = VCV_ES4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat4)
summary(mod_ES6) ##
## Multivariate Meta-Analysis Model (k = 87; method: REML)
##
## logLik Deviance AIC BIC AICc
## -41.2407 82.4815 92.4815 104.6947 93.2410
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0213 0.1459 27 no Study_ID
## sigma^2.2 0.0304 0.1744 87 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 85) = 293.4047, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 85) = 5.6275, p-val = 0.0051
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute -0.1285 0.1092 -1.1760 85 0.2429 -0.3457 0.0887
## Stress_durationChronic 0.1644 0.0523 3.1420 85 0.0023 0.0603 0.2684
##
## Stress_durationAcute
## Stress_durationChronic **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES6) ## R2_marginal R2_coditional
## 0.2226789 1.0000000# Orchard plot
orchard_plot(mod_ES6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Social enrichment
Does EE also include a manipulation of social environment (i.e., number of individuals in EE relative to control)?
VCV_ES5 <- make_VCV_matrix(data = dat5, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES7<- rma.mv(yi = lnRR_ESa, V = VCV_ES5, mod = ~EE_social-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_ES7)##
## Multivariate Meta-Analysis Model (k = 90; method: REML)
##
## logLik Deviance AIC BIC AICc
## -44.6612 89.3224 99.3224 111.7091 100.0541
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0237 0.1540 29 no Study_ID
## sigma^2.2 0.0283 0.1683 90 no ES_ID
## sigma^2.3 0.0385 0.1961 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 319.5169, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 88) = 1.0646, p-val = 0.3493
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_socialNon-social 0.1083 0.1226 0.8838 88 0.3792 -0.1353 0.3519
## EE_socialSocial 0.1791 0.1235 1.4501 88 0.1506 -0.0664 0.4246
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES7) ## R2_marginal R2_coditional
## 0.01313274 0.58063226# Orchard plot
orchard_plot(mod_ES7, mod = "EE_social", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Exercise enrichment
Does the form of enrichment include exercise through a running wheel or treadmill?
mod_ES8<- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES8)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -45.5487 91.0975 101.0975 113.5965 101.8118
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0337 0.1837 30 no Study_ID
## sigma^2.2 0.0273 0.1653 92 no ES_ID
## sigma^2.3 0.0049 0.0702 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 309.5719, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 1.9734, p-val = 0.1450
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1136 0.0762 1.4910 90 0.1395 -0.0378 0.2649
## EE_exerciseNo exercise 0.1670 0.1002 1.6672 90 0.0989 -0.0320 0.3660
##
## EE_exerciseExercise
## EE_exerciseNo exercise .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES8) ## R2_marginal R2_coditional
## 0.009818664 0.926129248# Orchard plot
orchard_plot(mod_ES8, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Age of enrichment
The age at which the individuals were exposed to environmental enrichment.
This needs to be redone after recategorising
mod_ES9 <- rma.mv(yi = lnRR_ESa, V = lnRRV_ES, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES9) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.9629 85.9259 97.9259 112.8577 98.9503
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0547 0.2340 30 no Study_ID
## sigma^2.2 0.0170 0.1305 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 299.7681, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.0451, p-val = 0.0329
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Age_EE_exposureAdult 0.1576 0.1030 1.5299 89 0.1296 -0.0471 0.3623
## Age_EE_exposureJuvenile -0.0146 0.1733 -0.0840 89 0.9332 -0.3589 0.3298
## Age_EE_exposureUnclear 0.1693 0.0650 2.6053 89 0.0108 0.0402 0.2985
##
## Age_EE_exposureAdult
## Age_EE_exposureJuvenile
## Age_EE_exposureUnclear *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES9) ## R2_marginal R2_coditional
## 0.03931516 0.99999999# Orchard plot
orchard_plot(mod_ES9, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
‘Pure effect sizes’
Comparing ES relative to stress we can see that the increase is about 18% but ES relative to control is 12% (i.e., this is the interaction).
We can also see that ES relative to control 9% so these individuals actually do better than baseline individuals. There is also no significant difference between ES and EC individuals (although estimate is slightly lower), meaning that enrichment in the presence of stress can result in learning and memory that is just as good as enriched individuals without stress
Enrichment relative to control
mod_E20 <- rma.mv(yi = lnRR_E2a, V = lnRRV_E2, random = list(~1|Study_ID,
# ~ 1|Strain, does not run as we have NA
~1|ES_ID),
test = "t",
data = dat)
summary(mod_E20) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -84.0627 168.1253 174.1253 181.6579 174.4012
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0119 0.1089 30 no Study_ID
## sigma^2.2 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 18.3889, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1272 0.0619 2.0546 91 0.0428 0.0042 0.2501 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1orchard_plot(mod_E20, mod = "Int", xlab = "lnRR")
Stress relative to control
mod_S20 <- rma.mv(yi = lnRR_S2a, V = lnRRV_S2, random = list(~1|Study_ID,
# ~ 1|Strain, does not run as we have NA
~1|ES_ID),
test = "t",
data = dat)
summary(mod_S20) ###
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -90.1569 180.3138 186.3138 193.8464 186.5897
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0167 0.1291 30 no Study_ID
## sigma^2.2 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 29.0591, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1112 0.0660 -1.6859 91 0.0952 -0.2422 0.0198 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1orchard_plot(mod_S20, mod = "Int", xlab = "lnRR")
Enrichment + stress relative to control
mod_ES20 <- rma.mv(yi = lnRR_ES2a, V = lnRRV_ES2, random = list(~1|Study_ID,
# ~ 1|Strain, does not run as we have NA
~1|ES_ID),
test = "t",
data = dat)
summary(mod_ES20) ###
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -88.6424 177.2848 183.2848 190.8174 183.5607
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0122 0.1103 30 no Study_ID
## sigma^2.2 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 22.5439, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.0934 0.0639 1.4621 91 0.1472 -0.0335 0.2204
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1orchard_plot(mod_ES20, mod = "Int", xlab = "lnRR")
Enrichment + stress relative to stress
mod_E30 <- rma.mv(yi = lnRR_E3a, V = lnRRV_E3, random = list(~1|Study_ID,
# ~ 1|Strain, does not run as we have NA
~1|ES_ID),
test = "t",
data = dat)
summary(mod_E30) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -91.7018 183.4037 189.4037 196.9362 189.6795
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0115 0.1072 30 no Study_ID
## sigma^2.2 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 26.3880, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1843 0.0629 2.9317 91 0.0043 0.0594 0.3092 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1orchard_plot(mod_E30, mod = "Int", xlab = "lnRR")
Enrichment + stress relative to enrichment
mod_S30 <- rma.mv(yi = lnRR_S3a, V = lnRRV_S3, random = list(~1|Study_ID,
# ~ 1|Strain, does not run as we have NA
~1|ES_ID),
test = "t",
data = dat)
summary(mod_S30) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -82.4100 164.8200 170.8200 178.3526 171.0958
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 30 no Study_ID
## sigma^2.2 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 10.9070, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0093 0.0503 -0.1855 91 0.8533 -0.1092 0.0906
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1orchard_plot(mod_S30, mod = "Int", xlab = "lnRR")
Combined orchard plot
#TODO need to figure out how to label y-axis. Also need to think about the most logical order to present these so that readers will understand what they are showing
mod_list <- list(mod_E20, mod_S20, mod_ES20, mod_E30, mod_S30)
mod_res <- lapply(mod_list, function(x) mod_results(x, mod = "Int"))
merged <- submerge(mod_res[[1]], mod_res[[2]], mod_res[[3]], mod_res[[4]], mod_res[[5]], mix = T)
#merged$mod_table$name2 <- c("EC vs CC", "CS vs CC", "ES vs CC", "ES vs EC", "ES vs CS")
orchard_plot(merged, mod = "Int", xlab = "lnRR", angle = 0) #+ geom_point(data = merged$mod_table, aes(y = name2, x = estimate))
Modelling with SMD
Stress
Intercept model
Why does the funnel plot look so strange and why is I2 higher?
mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0a) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -118.9545 237.9091 245.9091 255.9525 246.3742
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.6193 0.7869 30 no Study_ID
## sigma^2.2 0.4872 0.6980 92 no ES_ID
## sigma^2.3 0.0000 0.0001 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 27513.4051, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.3674 0.1666 -2.2044 91 0.0300 -0.6984 -0.0363 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_S0a) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.979242e-01 5.585262e-01 4.393980e-01 2.491240e-09funnel(mod_S0a) 
funnel(mod_S0a, yaxis="seinv",) 
Environmental Enrichment
Intercept model
Why does the funnel plot look so strange and why is I2 higher?
mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -116.8436 233.6871 241.6871 251.7305 242.1522
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3557 0.5964 30 no Study_ID
## sigma^2.2 0.5325 0.7297 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 14139.3073, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.7058 0.1390 5.0791 91 <.0001 0.4298 0.9818 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_E0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.974155e-01 3.994255e-01 5.979900e-01 1.844209e-09funnel(mod_E0a) 
Interaction
Intercept model
Why does the funnel plot look so strange and why is I2 higher?
mod_ES0a <- rma.mv(yi = SMD_ESa, V = VCV_ESa, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -128.8823 257.7646 265.7646 275.8080 266.2297
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.8719 0.9338 30 no Study_ID
## sigma^2.2 0.5852 0.7650 92 no ES_ID
## sigma^2.3 0.0000 0.0001 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 12257.6601, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.6295 0.1952 3.2255 91 0.0017 0.2418 1.0172 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_ES0a) #I2 is very high! 99.37%## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.937219e-01 5.946336e-01 3.990883e-01 6.436368e-09funnel(mod_ES0a)
Social enrichment
Does EE also include a manipulation of social environment? Note that we excluded any studies that exclusively used social enrichment.s